function Step7generate_costShocks_35X31()
% *************************************************************************
% GENERATING COST SHOCKS - Higher-order input linkages: 35 X 31: all sectors 31 countries
% *************************************************************************
%       (1) This code collects the following matrices:
%           - PPI
%           - WIOD related matrices: GAMMA, D, TOWEIGHTS
%           - Exchange rate
%       (2) Checks country, sector or country-sector alignments across matrices
%       (3) Generates special-sized GAMMA matrices 
%       (4) Generates cost shocks under
%           Cost shock: higher-order input linkages
%       (5) Aggregated PPI and cost shocks are generated
%       (6) Output is converted to a table & saved in
%           - country-sector level       : 'costShocks_allPeriods_table_higherOrder.xlsx'
%           - aggregated to country level: 'aggregatedCostShocks_allPeriods_table_higherOrder.xlsx'
% *************************************************************************
clear; clc;
tic;

%% Settings
rootfolder      = input('Enter the location for the rootfolder in single quotes (''''): \n');
projectfolder   = [rootfolder '\analysis'];
ppifolder       = [projectfolder '\PPI\MatlabFiles\'];
gammafolder     = [projectfolder '\WIOD\GAMMA\MatlabFiles\'];
alphafolder     = [projectfolder '\WIOD\ALPHA\MatlabFiles\'];
sigmafolder     = [projectfolder '\WIOD\SIGMA\MatlabFiles\'];
sweightsfolder  = [projectfolder '\WIOD\SWEIGHTS\MatlabFiles\'];
toweightsfolder = [projectfolder '\WIOD\TOWEIGHTS\MatlabFiles\'];
fxfolder        = [projectfolder '\FX\MatlabFiles\'];
outputfolder    = [projectfolder '\COSTSHOCKS\35X31\'];
betafolder      = [rootfolder '\data\6_BETA\'];
if exist(outputfolder,'dir') == 0
    mkdir(outputfolder);
end

cd(projectfolder);

% Other input files
beta_file = [betafolder 'sectorspecificbetas.xlsx'];
%% -- (Start looping over the year & month)
fprintf([' * Generating cost shocks (35 X 31) for the year: \n']);

for y = 1995:2011
    yr = num2str(y);
    fprintf([' - ' yr ' \n']);

    for m = 1:12
        mo = num2str(m);
        
%         fprintf([' ************************************** \n  * Running year: ' yr ', month: ' mo '  \n ************************************** \n']);
        
        %% Download vectors & matrices
        % one for big services scenarios & full linkages
        
        %         % Download Gamma_all matrix
        %         GAMMAdata = importdata([gammafolder,'WIOD_gamma_31cty35sec_' yr '.csv']);
        %         Gamma_text = GAMMAdata.textdata(2:end,1:3);
        %         GAMMA = GAMMAdata.data;
        
        % Download Gamma_OO matrix
        GAMMAoodata = importdata([gammafolder,'WIOD_gamma_31cty35sec_' yr '_OO.csv']);
        Gammaoo_text = GAMMAoodata.textdata(2:end,1:3);
        GAMMAoo = GAMMAoodata.data(:,1:end-2);
        
        % Download Gamma_ON matrix
        GAMMAondata = importdata([gammafolder,'WIOD_gamma_31cty35sec_' yr '_ON.csv']);
        Gammaon_text = GAMMAondata.textdata(2:end,1:3);
        GAMMAon = GAMMAondata.data;
        
        % Download Gamma_NO matrix
        GAMMAnodata = importdata([gammafolder,'WIOD_gamma_31cty35sec_' yr '_NO.csv']);
        Gammano_text = GAMMAnodata.textdata(2:end,1:3);
        GAMMAno = GAMMAnodata.data;
        
        % Download Gamma_NN matrix
        GAMMAnndata = importdata([gammafolder,'WIOD_gamma_31cty35sec_' yr '_NN.csv']);
        Gammann_text = GAMMAnndata.textdata(2:end,1:3);
        GAMMAnn = GAMMAnndata.data;
        
        % List of countries
        listoo_country = unique(Gammaoo_text(:,2));
        liston_country = unique(Gammaon_text(:,2));
        listno_country = unique(Gammano_text(:,2));
        listnn_country = unique(Gammann_text(:,2));
        assert(sum(strcmp(listoo_country,liston_country)) == size(listoo_country,1));
        assert(sum(strcmp(liston_country,liston_country)) == size(listno_country,1));
        assert(sum(strcmp(listoo_country,listno_country)) == size(listnn_country,1));
        
        list_country = listoo_country;
        n_country = size(list_country,1);
        
        % List of sectors - CAREFUL! about the order of sectors
        list_sector = [unique(Gammaoo_text(:,3)); unique(Gammano_text(:,3))];
        list_sector2 = [unique(Gammaon_text(:,3)); unique(Gammann_text(:,3))];
        assert(sum(strcmp(list_sector,list_sector2)) == size(list_sector,1));
        n_sector = size(list_sector,1);
        n_cs = n_country * n_sector;
        
        % Countries & sectors for OO
        list_country_oo = unique(Gammaoo_text(:,2));
        n_country_oo = size(list_country_oo,1);
        list_sector_oo = unique(Gammaoo_text(:,3));
        n_sector_oo = size(list_sector_oo,1);
        n_cs_oo = n_country_oo * n_sector_oo;
        
        % Countries & sectors for NN
        list_country_nn = unique(Gammann_text(:,2));
        n_country_nn = size(list_country_nn,1);
        list_sector_nn = unique(Gammann_text(:,3));
        n_sector_nn = size(list_sector_nn,1);
        n_cs_nn = n_country_nn * n_sector_nn;
        
        % Generate identity & 1 & 0 matrices
        Ioo = eye(size(GAMMAoo));
        Ion = eye(size(GAMMAon));
        Ino = eye(size(GAMMAno));
        Inn = eye(size(GAMMAnn));
        
        % Generate D matrix - for oo
        Doo = diag(1 - sum(GAMMAoo,2));
        Doo_inv = Doo \ Ioo;
        
        % Download PPI vector
        ppioodata = importdata([ppifolder,'PPI_' yr mo '.csv']);
        ppioo_text = ppioodata.textdata(2:end,1:5);
        PPIoo = ppioodata.data;
        % -size check
        assert(size(PPIoo,1) == n_cs_oo )
        assert(size(PPIoo,2) == 1 )
        
        % Download E matrix versus USD - short version
        Edata = importdata([fxfolder,'FX_' yr mo '.csv']);
        E_text = Edata.textdata(2:end,1:4);
        E_short = Edata.data;
        % -size check
        assert(size(E_short,1) == n_country)
        assert(size(E_short,2) == 1 )
        
        % Compute E_0 matrix - Kronecker product: E*1's
        E_0 = kron(E_short,ones(n_sector,1));
        
        % Download BILATERAL E matrix - short version
        E_bil_data = importdata([fxfolder,'FX_bil_' yr mo '.csv']);
        E_bil_text = E_bil_data.textdata(2:end,1:4);
        E_bil_short = E_bil_data.data;
        % -size check
        assert(size(E_bil_short,1) == n_country^2)
        assert(size(E_bil_short,2) == 1 )
        
        % Compute E_bil matrix - Kronecker product: E*1's
        E_bil_o = kron(E_bil_short,ones(n_sector_oo,1));
        E_bil_n = kron(E_bil_short,ones(n_sector_nn,1));
        E_bil = [E_bil_o; E_bil_n];
        E_bil2 = kron(E_bil_short,ones(n_sector,1));
        
        % Download total output weights
        TOWEIGHTSdata = importdata([toweightsfolder,'WIOD_TOweights_' yr '.csv']);
        TOweights_text = TOWEIGHTSdata.textdata(2:end,1:3);
        TOWEIGHTS = TOWEIGHTSdata.data;
        
%         fprintf('  --- 35 X 31 matrices downloaded (18X31 and 17X31 as well) \n');
        
        %% Generate special-sized GAMMAs
        % Generate BIG GAMMA
        GAMMA = [GAMMAoo GAMMAon; GAMMAno GAMMAnn];
        Gamma_text = [Gammaoo_text; Gammano_text];
       
        gammaoo_country = Gammaoo_text(:,2);
        gammann_country = Gammann_text(:,2);
        
        GAMMAo = [GAMMAoo GAMMAon];
        gammao_text = Gammaoo_text;
        gammao_country = Gammaoo_text(:,2);
        gammao_sector = Gammaoo_text(:,3);

        GAMMAn = [GAMMAno GAMMAnn];
        gamman_text = Gammano_text;
        gamman_country = Gammann_text(:,2);
        gamman_sector = Gammann_text(:,3);
        
        % Generate D
        D = diag(1 - sum(GAMMA,2));
        
        % Generate special-sized GAMMAs to multiply with special-sized E_bil
        % - GAMMA_oo_tilde
        list = '';
        for i = 1:n_country
            eval(['prt_' num2str(i) ' = GAMMAoo((i-1)*n_sector_oo+1 : i*n_sector_oo,:);']) ;
            if i == n_country
                list = [list 'prt_' num2str(i)];
            else
                list = [list 'prt_' num2str(i) ','];
            end
        end
        
        GAMMAoo_tilde = blkdiag(prt_1,prt_2,prt_3,prt_4,prt_5,prt_6,prt_7,prt_8,prt_9,prt_10,prt_11,prt_12,prt_13,prt_14,prt_15,prt_16,prt_17,prt_18,prt_19,prt_20,prt_21,prt_22,prt_23,prt_24,prt_25,prt_26,prt_27,prt_28,prt_29,prt_30,prt_31);
        clear prt*
        gammaoo_tilde_country_horizontal = repmat(gammaoo_country,n_country,1);
        
        % - GAMMA_on_tilde
        list = '';
        for i = 1:n_country
            eval(['prt_' num2str(i) ' = GAMMAon((i-1)*n_sector_oo+1 : i*n_sector_oo,:);']) ;
            if i == n_country
                list = [list 'prt_' num2str(i)];
            else
                list = [list 'prt_' num2str(i) ','];
            end
        end
        
        GAMMAon_tilde = blkdiag(prt_1,prt_2,prt_3,prt_4,prt_5,prt_6,prt_7,prt_8,prt_9,prt_10,prt_11,prt_12,prt_13,prt_14,prt_15,prt_16,prt_17,prt_18,prt_19,prt_20,prt_21,prt_22,prt_23,prt_24,prt_25,prt_26,prt_27,prt_28,prt_29,prt_30,prt_31);
        clear prt*
        gammaon_tilde_country_horizontal = repmat(gammann_country,n_country,1);
        
        % * GAMMA_o_tilde
        GAMMAo_tilde = [GAMMAoo_tilde GAMMAon_tilde];
        
        % - GAMMA_no_tilde
        list = '';
        for i = 1:n_country
            eval(['prt_' num2str(i) ' = GAMMAno((i-1)*n_sector_nn+1 : i*n_sector_nn,:);']) ;
            if i == n_country
                list = [list 'prt_' num2str(i)];
            else
                list = [list 'prt_' num2str(i) ','];
            end
        end
        
        GAMMAno_tilde = blkdiag(prt_1,prt_2,prt_3,prt_4,prt_5,prt_6,prt_7,prt_8,prt_9,prt_10,prt_11,prt_12,prt_13,prt_14,prt_15,prt_16,prt_17,prt_18,prt_19,prt_20,prt_21,prt_22,prt_23,prt_24,prt_25,prt_26,prt_27,prt_28,prt_29,prt_30,prt_31);
        clear prt*
        gammano_tilde_country_horizontal = repmat(gammaoo_country,n_country,1);

        
        % - GAMMA_on_tilde
        list = '';
        for i = 1:n_country
            eval(['prt_' num2str(i) ' = GAMMAnn((i-1)*n_sector_nn+1 : i*n_sector_nn,:);']) ;
            if i == n_country
                list = [list 'prt_' num2str(i)];
            else
                list = [list 'prt_' num2str(i) ','];
            end
        end
        
        GAMMAnn_tilde = blkdiag(prt_1,prt_2,prt_3,prt_4,prt_5,prt_6,prt_7,prt_8,prt_9,prt_10,prt_11,prt_12,prt_13,prt_14,prt_15,prt_16,prt_17,prt_18,prt_19,prt_20,prt_21,prt_22,prt_23,prt_24,prt_25,prt_26,prt_27,prt_28,prt_29,prt_30,prt_31);
        clear prt*
        gammann_tilde_country_horizontal = repmat(gammann_country,n_country,1);
        
        % * GAMMA_n_tilde
        GAMMAn_tilde = [GAMMAno_tilde GAMMAnn_tilde];
        
        % * GAMMA_tilde
        GAMMA_tilde = [GAMMAo_tilde; GAMMAn_tilde];
        
        Gamma_tilde_country_horizontal = [gammaoo_tilde_country_horizontal; gammaon_tilde_country_horizontal];
        Gamma_tilde_country_horizontal2 = [gammano_tilde_country_horizontal; gammann_tilde_country_horizontal];
        
        assert(sum(strcmp(Gamma_tilde_country_horizontal,Gamma_tilde_country_horizontal2)) == size(Gamma_tilde_country_horizontal,1));
        
        %% Check country - sector alignment across matrices
        
        % Sector list: Gamma versus PPI
        gammaoo_sector = Gammaoo_text(:,3);
        ppioo_sector   = ppioo_text(:,4);
        assert(sum(strcmp(ppioo_sector,gammaoo_sector)) == size(gammaoo_sector,1));
        
        % Sector list: total output weights versus Gamma
        toweights_sector = TOweights_text(:,3);
        assert(sum(strcmp(toweights_sector,gammaoo_sector)) == size(gammaoo_sector,1) );
        
%         fprintf('  --- Sector list: aligned - 35 X 31 \n');
        
        % Country list: Gamma versus PPI
        gammaoo_country = Gammaoo_text(:,2);
        ppioo_country = ppioo_text(:,5);
        assert(sum(strcmp(ppioo_country,gammaoo_country)) == size(gammaoo_country,1));
        
        % Country list: E_bil versus Gamma_tilde
        E_bil_short_country = E_bil_text(:,4);
        E_bil_country = [repelem(E_bil_short_country,n_sector_oo,1); repelem(E_bil_short_country,n_sector_nn,1)];
        assert(sum(strcmp(E_bil_country,Gamma_tilde_country_horizontal)) == size(Gamma_tilde_country_horizontal,1));
        
        % Country list: total output weights versus Gamma
        toweights_country = TOweights_text(:,2);
        assert(sum(strcmp(toweights_country,gammaoo_country)) == size(gammaoo_country,1) );
        
%         fprintf('  --- Country list: aligned - 35 X 31 \n');

        
        %% --- Generate cost shocks ---
        
        %% Cost shock: for higher-order input linkages
        
        % Compute cost shock
        IminusGAMMAnn_inv = (Inn - GAMMAnn) \ Inn;
        costho_firstTerm = Doo_inv * (Ioo - GAMMAoo -  GAMMAon * IminusGAMMAnn_inv * GAMMAno ) * PPIoo ;
        costho_secondTerm =  Doo_inv * (GAMMAo_tilde + (GAMMAon * IminusGAMMAnn_inv * GAMMAn_tilde) ) * E_bil;
        COST_ho = costho_firstTerm - costho_secondTerm;
        
%         fprintf('  --- Generated Cost shocks for higher-order input linkages - for 35 X 31 \n');
        
        %% Gather PPI & cost shocks in a table
        headers_textdata = ppioodata.textdata(1,1:5);
        textdata = ppioodata.textdata(2:end,1:5);
        
        % Get the list of all generated cost shocks, and add PPI
        if y == 1995 && m == 1
            numdatanames = who('-regexp','PPI*');
            numdatanames = [numdatanames; who('-regexp','COST*')];
        end
        
        numdata = []; headers_numdata = '';
        for n = 1:length(numdatanames)
            numdata = [numdata eval(numdatanames{n})];
            headers_numdata = [headers_numdata cellstr(numdatanames{n})];
        end
        
        %headers_numdata = ['PPI' 'COST1' ];
        %numdata = [PPI COST1 ];
        
        TABLE = [];
        for i = 1: length(headers_textdata)
            clear T;
            T = table(textdata(:,i),'variableNames',cellstr(headers_textdata{i}));
            TABLE = [TABLE T];
        end
        
        for i = 1: length(headers_numdata)
            clear T;
            T = table(numdata(:,i),'variableNames',cellstr(headers_numdata{i}));
            TABLE = [TABLE T];
        end
        
        %% Generate aggregated PPI & cost shocks at the country level
        % - using 2002 total output weights
        
        % Download total output weights 2002
        TOWEIGHTSdata2002 = importdata([toweightsfolder,'WIOD_TOweights_2002.csv']);
        TOweights2002_text = TOWEIGHTSdata2002.textdata(2:end,1:3);
        TOWEIGHTS2002 = TOWEIGHTSdata2002.data;
        TOWEIGHTS2002_tomultiply = reshape(TOWEIGHTS2002,[n_sector_oo,n_country_oo]);
        
        % Reshape for multiplying by weighting matrix
        for n = 1:length(numdatanames)
            eval([ numdatanames{n} '_forweighing = transpose(reshape(' numdatanames{n} ',[n_sector_oo,n_country_oo]));']);
        end
        
        % Multiply and take the diagonal: to weigh
        for n = 1:length(numdatanames)
            eval([ numdatanames{n} '_ag = diag(' numdatanames{n} '_forweighing * TOWEIGHTS2002_tomultiply);']);
        end
       
        
        %% Gather aggregated PPI & cost shocks in a table
        % Get the list of generated aggregated cost shocks, and add PPI
        numdata_aggregate = []; headers_numdata_aggregate = ''; % initialize
        for n = 1:length(numdatanames)
            numdata_aggregate = [numdata_aggregate eval([numdatanames{n} '_ag'])];
            headers_numdata_aggregate = [headers_numdata_aggregate cellstr([numdatanames{n} '_ag'])];
        end
        
        %headers_numdata_aggregate = list2Cell('PPI_aggregate,COST1_aggregate);
        %numdata_aggregate = [PPI_aggregate COST1_aggregate];
        
        % Prepare text data
        selectC = [1:n_sector_oo:n_cs_oo];
        country_text = Gammaoo_text(selectC,2);
        year_text = cellstr(repmat(yr,size(country_text)));
        month_text = cellstr(repmat(mo,size(country_text)));
        
        % Generate table
        Tcountry = table(country_text,'variableNames',{'country'});
        Tyear    = table(year_text,'variableNames',{'year'});
        Tmonth   = table(month_text,'variableNames',{'month'});
        AGGTABLE = [Tcountry Tyear Tmonth];
        
        for i = 1: length(headers_numdata_aggregate)
            clear T;
            T = table(numdata_aggregate(:,i),'variableNames',cellstr(headers_numdata_aggregate{i}));
            AGGTABLE = [AGGTABLE T];
        end
        
        %% Stack tables for each period in a large table with all periods
        if y == 1995 && m == 1
            TABLE_allP = []; AGGTABLE_allP = []; % initialize
        end
        TABLE_allP = [TABLE_allP; TABLE]; % cs breakdown
        AGGTABLE_allP = [AGGTABLE_allP; AGGTABLE]; % c breakdown
        
        %% -- (End the loop for the month)
    end
    %% -- (End the loop for the year)
end

%% Export to excel/csv
% Cost shock for each country-sector
output_filename = ['costShocks_allPeriods_table_higherOrder.xlsx'];
writetable(TABLE_allP,fullfile(outputfolder,output_filename),'WriteRowNames',true);
fprintf(['  --- Cost shocks exported to excel: ' output_filename '  \n']);

% Aggregated cost shock for each country
output_filename = ['aggregatedCostShocks_allPeriods_table_higherOrder.xlsx']; % The file used in factor analysis
writetable(AGGTABLE_allP,fullfile(outputfolder,output_filename),'WriteRowNames',true);
fprintf(['  --- Aggregated cost shocks exported to excel: ' output_filename '  \n']);

toc
pause(3)
exit
end
